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P,: ABSTRACT 
(D 

' PSR B1937+21 was the first millisecond pulsar ever measured. Which has 

^ ■ been appeared as a singular pulsar. The high-precision observation of this pulsar 

■ shows systemic long-term variation in the residual of time of arrivals. This paper 

^ • modelled the secular variation by the orbital precession induced time delay of a 

O ■ binary pulsar system. The fitting requires that as a binary pulsar, PSR B1937+21 

Q\ • should have small companion star, m2 ~ 10~'^Mq, and projected semi-major axis, 
X ~ 10~^s to lO^^s. Which corresponds to ignorable radial velocity of the pulsar 



O 



O ■ to the line of sight. This might explain why it has been measured as a singular 

^ ■ pulsar instead of a binary pulsar. 

2 ' Subject headings: gravitation — pulsar: individual (PSR B1937+21) — relativity- 

^ ■ stars: binaries-stars: fundamental parameters — stars: neutron 

> 

X 

;h ' 1. Introduction 

PSR B1937+21, the first millisecond pulsar and also the fastest rotating pulsar up to 
date, was discovered by Backer & Kulkarni et al. (1982). Which has been measured as a 
singular pulsar. The systemic long-term variation in the residual of time of arrivals (TOAs) 
has been measured by different authors, Kaspi et al. (1994), Lommen & Backer (2001). 
The variation may be caused by planet around the pulsar. However, such possibility was 
constrained stringently by Thorsett & Phillips (1992). 

This paper provides an alternative model to interpret the systematic variation in TOAs, 
the orbital precession induced time delay of a binary pulsar system. Three numerical solu- 
tions are obtained through fitting the residual of TOAs, all of which demands small compan- 
ion star mass, ~ lO^^M©, and in turn small projected semi-major axis, x ~ 10~^ or 10^^. 
Therefore, the corresponding radial velocity of the pulsar to the line of sight (LOS) is so 
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small that it cannot be resolved from the timing noise of TOAs. Which might explain why 
the companion star has not been observed from timing measurement. 

This paper shows that although the companion star of PSR B 1937+21 cannot be mea- 
sured directly from the velocity curve, however, the orbital parameters could still be extracted 
through the long-time variation of TOAs. 

By the orbital precession model, when the fitted orbital parameters are put into the 
standard timing model, the behavior of PSR B1937+21 as binary pulsar should be similar to 
PSR J205 1-0827 (Doroshenko et al. 2001) and PSR B1957+21 (Arzoumanian et al. 1994), 
which have similar companion masses and orbital period. Moreover, when the additional 
time delay is also added to the standard timing model, the systemic variation of TOAs of 
PSR B1937+21 should be eliminated, so are PSR J2051-0827 and PSR B1957+21. PSR 
B1937+21 might be another millisecond pulsar whose companion star is evaporated by the 
radiation of the pulsar as PSR B1957+21. 

2. The orbital precession model 

Barker & O'Connell (1975) gives the first gravitational two-body equation with spins, 
in which the orbital angular momentum, vector, L, is expressed as precesses around a vector 
combined by Si and S2, and the magnitude is small 2 Post Newtonian (PN), typically 
1° • 10^''yr~^, which is insignificant in pulsar timing measurement. 

On the other hand, the precession of orbital angular momentum vector, L, around the 
total angular momentum vector, J, in the special case that ignores the spin angular momen- 
tum of one star in a binary pulsar system, has been studied by many authors (Apostolatos 
et al. 1994, Kidder 1995, Wex & Kopeikin 1999). In this case the precession velocity of 
the orbit plane is significant, 1.5 PN, typically l°yr~^ for binary pulsars with orbital period 
around 10 hours. Which is measurable in pulsar timing. 

Therefore, the two kinds of expressions on the precession velocity of L seems contrary 
each other. To investigate their relationship, the first and most important question to be 
answered is: which vector should L precess around ? The two kinds of expressions both agree 
that Si and S2, as well as the vector combined by them precess rapidly (1.5PN) around the 
total angular momentum, J. Since J is at rest relative to LOS (after counting out the 
proper motion), so the small precession velocity of L around the vector combined by Si and 
S2 doesn't mean that precession velocity of L around J is also small, and in turn relative 
velocity LOS is small. 
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The situation is analogy to the following case. If the binary system we observe is 
replaced by a solar system, one cannot say that the velocity of planet A is ignorable because 
its velocity with respect to planet B is very small. Instead only when the velocity of planet A 
is very small relative to the Baryon center of the solar system which is at rest to the observer 
is very small, can one conclude that this velocity is ignorable. 

Obviously the role of J in a binary system is equivalent that baryon center in the solar 
system. And only the velocity relative to J make sense in observation. Therefore, the 
precession of L should be expressed as around J to compare with observation. 

In the one spin case the precession velocity of L must be around J and the magnitude 
must be significant 1.5PN (Apostolatos et al. 1994, Kidder 1995, Wex & Kopeikin 1999). 
Gong develop the one spin to the general two spins case, and obtains that still the precession 
velocity of L must be around J and must be significant 1.5PN (Gongl 2003). 

The difference is that in the two spins case, the magnitude of the precession velocity of 
Oq varies rapidly due to the variation of the sum the spin angular momenta of the two stars, 
S, which can lead to significant secular variabilities in binary pulsars. Whereas, the one spin 
case predicts a constant magnitude of fio, and thus it cannot explain the significant secular 
variabilities measured in binary pulsars. 

Since the spin angular momenta of the two stars in binary pulsar, 5*1 and 5*2, are much 
smaller than that of the orbital angular momentum, L, therefore, the sum of the two spins, 
S = Si + 82, is also much smaller than L, or, S/ L must be satisfied for a general binary 
pulsar system. Moreover since J = L + S (J, L and S forms a triangle), the constraint, 
S/L <^ 1, means the misalignment angle between J and L must be very small, Xlj <^ 1. 
The conservation of the total angular momentum of a binary pulsar system, J = 0, leads to 
(Barker & O'Connell 1975), 

Oq X L = — Oi X Si — 02 X S2 , (1) 

the magnitude of the left hand side of Eq(l) is Jlo-Lsin A^j. Prom which one can see that for 
a given torque caused by and 5*2 at the right hand side of Eq(l), there are two choices 
for Qq. The first one is that Qq is small in the case that Xlj ~ 1, the second one is that Qq 
is large in the case that Xlj <^ 1. As analyzed above, not matter from the point of view of 
observation, or from the constraint that must be satisfied by a binary pulsar {S/L <^ 1), Qq 
have to follow the second choice. Which means L must preccss around J with a very small 
angle of precession cone, Alj <^ 1, and the magnitude of flo must be significant. 

By Eq(l), the velocity of the orbital precession, Qq, which is along J and including two 
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spins can be obtained (Gongl 2003). 

9,0 = p9,2 sin Xls + pi^i - Q2) ^ sin A^^^ , (2) 

where p — ^.^^^^ , and — Si cos 77551 denotes the component of Si in the plane determined 

by S and J. X\g_^ is the angle between sf and L. Note that LsinA^j ~ pS is used in Eq(2), 
since S/L <^ 1. The right-hand side of Eq(2) can as well be written by replacing subscribes 
1 with 2 and 2 with 1. The precession velocity of the two stars, Jli and Q2, are given by 
Barker & O'Connell's equation. 

Actually Eq(2) can be regarded as transforming the direction of Qq of the Barker & 
O'Connell's equation (two spins) to J; and it can also be regarded as add one more spin into 
the former one spin case given by Apostolatos et al. (1994), Kidder (1995), Wex & Kopeikin 
(1999). 

The scenario of the motion of a binary system has been discussed by authors, Smarr & 
Blandford (1976), Hamilton & Sarazin (1982). In which L, Si and S2 all precess around J 
rapidly (1.5PN). And the relative velocities of L to Si and S2 are very small. From which a 
slow precession of L around J is impossible because of the triangle constraint, J = L + S. 
The scenario can solve all the puzzle concerning the motion of the three vectors, L, Si and 
S2. 

Since Si and S2 process with different velocities, Qi and Q2 respectively (mi 7^ 1TI2), 
then S varies in both magnitude and direction (Si, S2 and S form a triangle), then from the 
triangle of S, L and J, in react to the variation of S, L must vary in direction (|L| =const), 
which means the variation of A^j (J is invariable). 

The change of Xlj means that the orbital plane tilts back and forth. In turn, both 

Al5 and Aj5 vary with time. Therefore, from Eq(2), the derivative of the rate of orbital 
precession can be given by, 

Qo^ pfl + pfl , (3) 

where fl = Qq/p and Q is given by 

Q = fi2^i2^3^4 ~ ^12^1(^01^2 + ^12^3) + ^12^1^5 ) (4) 

where = ^i - (^2, ^^oi = ^^i - ^o, ^1 = ^sinA^, X2 = tanr?,,i, X3 = '""l^'Z^l^ , 
X4 = £2^!^ X5 = cot a" a" and p = - ^.^ff In which a = sin Xjs + £2£^, 
Svi = Si sin Aj5i and Sv2 = 'S'2 sin Aj52 represent components of Si and S2 that are vertical 
to J. 
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Note that Qi and Q2 are unchanged when ignoring the orbital decay (2.5PN), and XlSc 
arc unchanged also, since they decay much slower than the orbital decay (Apostolatos et al. 
1994). Through Eq(3), can be easily obtained (Gongl 2003) . 

The S-L couphng in a binary pulsar system results the following effects. (1) the apsidal 
motion of the orbital plane which can be absorbed by u; (Smarr & Blandford 1976) 

^obs ^ ^GR ^ ^ ^GR ^ _ ^5^) 

Notice that is a function of time due to the variation of S, as shown in Eq(2), whereas, 
uj'^^, the GR prediction of the advance of periastron is a constant. 

(2) the precession of the orbital plane which can be absorbed by x, 

X = —xVLq sin A lj sin r/o cot i . (6) 

of Eq(6) is usually much larger than which is caused by the gravitational wave 

induced orbital decay. 

(3) the nutation of the orbital plane which can be absorbed by P},. The variation in 
the precession velocity of the orbit results in a variation of orbital frequency (f^ = 2t:/ Pj^), 
1/'^ — Uh — 0,0 At . Then z>6 = Qq, which lead to the variation of Pj, (Gongl 2003), 

h^-^. (T) 

Notice that the coupling of the spin induced quadrupole moment with the orbit (Q-L cou- 
pling) can also cause apsidal motion which can explain the secular variabilities on cu and 

X measured in binary pulsars. However, such effect cannot explain Pf„ and the second and 

(3) 

third order of derivatives of parameters like, x, Pj, and P^; . Since the Q-L coupling is similar 
to that of the S-L coupling in the one spin case, which corresponds to a static precession or 
apsidal motion of the orbital plane. In other words, the velocity, Qq is a constant in these 
two cases. Whereas, in the two spins S-L coupling, Qq is a function of time, as given by 
Eq(2), which can not only explain the parameters that can be explained by the Q-L coupling 
or the one spin S-L coupling, but also explain parameters that cannot be explained by them. 



3. Application to PSR B1937-|-21 

The essential transformation relating solar system baryccntric time tb to pulsar proper 
time T is given by the equation (Manchester and Taylor 1977) 

h-to^T + An + AE + As + AA, (8) 
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where Ar is the " Roemer time delay" , which corresponds to the propagation time across the 
binary orbit; and A5 arc the orbital Einstein and Shapiro delays; and is a time delay 
related with aberration caused by rotation of the pulsar. The dominant term concerning the 
orbital effects is the Roemer delay, Ar, given by (Taylor & Weisberg 1989) 

Ar = xF{u; + u) , (9) 

where 

F(u + u) = sina;[cosM — e(l + 6r)] + [1 — e^(l — Sg)'^]^^'^ cosa;sin-u , (10) 

where u is the eccentric anomaly, and uj is the angular distance periastron from the node. 
In calculation, the small quantities dr and Sq due to aberration are ignored. 

When the contribution of orbital precession to the time of arrival is considered, the 
dominant term that absorbs the additional time delay is also the Roemer delay. Which can 
deviate from the standard one, A^^, and therefore, leads to an additional time delay, 

A = xF{cu' + u') - xF{cu + u) = 2x sin((]ot/2) cos ^ xflot cos 0. (11) 

where 9 = (uj + QQ)t/2 + uji. The additional time delay, A, given by Eq(ll), is a function 
of time due to the variation of Qq. Which can contribute the TOAs via Eq(8). And such 
dynamic effect might be interpreted as effect of propagation, like DM in some binary pulsars. 
The effect of l^o, ^0 can be absorbed by such parameters (3), DM, DM^^\ 

p(3) p 

etc (Gong2, Gong3, 2003). The orbital precession model predicts that -i— ~ ~ f22 for a 
same binary pulsar, so is the derivatives of x and DM. Evidences of such correlation can 
be found in binary pulsars, such as PSR J2051-0827 (Doroshenko et al. 2001) and PSR 
B1957+21 (Arzoumanian et al. 1994), PSR B1534+12 (Konacki et al. 2003). Further more 
the orbital precession model predicts that binary pulsars with long orbital period, i.e., for 
days, correspond to much smaller VLq and and in turn much smaller derivatives of x, 
and DM by Eq(6), Eq(7) and Eq(ll). Which make the corresponding derivatives of such 
binary pulsar more difficult to observe than that of binary pulsars with small P;,, i.e., a few 
hours. Such correlation can also be find in binary pulsars 

As shown in Fig 1, the amplitude of the measured time delay is about 10~®s. By simple 
estimation, if xVLqcos0 ~ 10^^''' in Eq(ll), then the orbital precession induced time delay 
can be A ^ lO^^s also. There are two possibilities to satisfy this relation: (1) small VLq and 
UJ like PSR B1855+09 (large P;,), and large x, however in such case, it will take much longer 
time to finish a variation hke that of Fig 1. And (2) large VLq and a;, i.e., ~ 10°yr~^, and 
small X, obviously case (2) can not only satisfy the amplitude but also the periodicity of 
variation. Therefore, the orbital precession model implies that PSR B1937+21 should have 
small X, and small (which corresponds to large ^-nd uj). 
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The vectors Si, S2 and S are studied in the coordinate system of the total angular 
momentum, in which the z-axis directs to J, and the x- and y-axes are in the invariance 
plane. S can be represented by Sp and Sy, the components parallel and vertical to the 
z-axis, respectively: 

S={Sv + Sp)'/^ . (12) 
Sp and Sv can be expressed (recall Syi, Sv2 and Sy form a triangle) as 

Sp = 5"! cos Xjsi + S2 cos XjS2 , 
Sv = (^^1 + ^^2 - 2SviSv2 COS775152)'/' , (13) 

where 775152 is the misalignment angle between Syi and Sv2, which can be written as 

VS1S2 = {^1 - ^2)t + . (14) 

Therefore, by the variation of S as function of time (in the case of one spin, S =const), one 
can obtain Qq as function of time through Eq(2). Thus the measured A°^* of Fig 1 can be 
fitted step by step via the orbital precession induced A given by Eq(ll), 

A'"'\h+i) - A^'^tk) = A{h+i) - A{h) . (15) 

Note that both Qq and cos^ are functions of time in Eq(ll). Which are responsible for 
variation of residual in Fig 1. 

Three numerical solutions can be obtained by the Monte Carlo solution of Eq(15). The 
solutions are shown in Table 1,3,5, respectively. The three solutions shows that Pb is from 
0.0708d to O.lOd, and x from 3.73-10-^s to 4.25-10-3s, and m2 from 0.0126Mq to O.O373M0. 
The possibility of other solutions cannot be excluded. 

The standard DD timing model uses 12 parameter concerning the orbit to fit the delays 
of Eq(8) (Taylor & Weisberg 1989). The total number of binary parameters of the orbital 
precession model as shown in Table 1 is also 12. There are some parameters in common, 
such as, e, Pb, uq and Tq. The difference is that the latter can cause an additional time delay, 
as shown in Eq(ll), through the parameters of the second row of Table 1 (including mi and 
1712 of the first row). Whereas, the former doesn't have this delay. In other words, the latter 
only adds one more delay term to the right hand side of Eq(8), and all other terms (effects) 
are the same as the former models. The relation of the orbital precession model with the 
DDGR model can be given by P^'" = P^^ + P^^, and by Eq(5). 

By the results of Table 1,3,5 and the orbital precession model, one can easily obtain 
the following Post-Keplerian parameters. Notice that Table 2,4,6 correspond to the solution 
1,2,3, respectively. 
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4. discussion 

The orbital precession model predicts that PSR B1937+21 is a binary pulsar with very 
small orbital period, Pi, = 1.7hr. Then there is a question, why it has been measured as a 
singular pulsar? Actually it can be answered by the fitted orbital parameters of Table 1. 
The radial velocity (to line of sight) of the pulsar is given (Shapiro & Teukolsky 1983), 

27rai sin i , , . . 
vi-n=-^^^^— -^yyy^[cos(u; + 'u) + ecosu;J , (16) 

where ai = ^^^^ . Then the observed pulsar period deviates from the intrinsic period by a 
small amount, 

C Pft 

Therefore, the contribution of the orbital precession effect to P°^* is about 5P^^ ~ 27rxP/Pb ~ 
1.6 • 10~^°s (corresponds to solution 1). Whereas, the last decimal of pulsar period, P, is 
4(2) • 10~^^s (Kaspi et al. 1994), which is approximately two order of magnitude larger than 
§pOP jj^ other words, the modulation of the radial velocity of the pulsar (induced by binary 
motion) to the TOAs is ignorable. Moreover the small orbital inclination angle {i ~ 10~^) 
given by solution 1 and 2, means that the orbital plane of the binary may be face on, then 
it is also impossible to observe the eclipse of this binary pulsar. These factor may make it 
difficult to observe the effects of the companion star. 

By the orbital precession model, a binary pulsar with very small orbital period, i.e., a 
few hours, can produce a long-term time delay i.e., years, as given by Eq(ll). Therefore, the 
companion star of PSR B1937+21, although unmeasurable in the orbital period time scale, 
due to small x, might imprint in the long-term measurement of TOAs. 

The results of the three numerical solutions can be tested by TEMPO. Put the orbital 
parameters, such as P^, x, e of Table 1 (3,5), as well as Post-Kepler parameters P;,, cj of Table 
2 (4,6) into the corresponding terms in DD model. Then behavior of PSR B1937+21 should 
be very similar to PSR B1957+20 and PSR J2051-0827, since they have similar companion 
mass and close orbital period. Moreover, when the orbital precession induced time delay of 
Eq(ll) is added into Eq(8), then the systematic variation of residuals of TOAs should be 
eliminated. Since different initial time corresponds to different phase angles, so the initial 
phases, hke Ui, (f)i may need to be adjusted when fitted by TEMPO. 

Assuming PSR J2051-0827 and PSR B1957-I-20 are fitted as a singular pulsar, then they 
should have similar the residual of TOAs (after counting out DM variation) as that of PSR 
B1937-I-21. And binary pulsars such as, PSR J0437-4715, with much larger P^, should take 
much longer time to finish the similar variation in residuals of TOAs hke PSR B1937-I-21, 
in the case that PSR J0437-4715 is fitted as a singular pulsar. 
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The evolutionary links between millisecond pulsars and their binary progenitor systems 
involve a number of exotic astrophysical phenomena. In their late stages of evolution, neutron 
stars in low-mass X-ray and pulsar binaries may evaporate their companions through the 
strength of their radiation (Alpar et al. 1982, Bhattacharya & van den Heuvel 1992, King 
et al 2003). PSR B 1957+20 has provided strong evidences to that scenario, the presumed 
binary pulsar PSR B1937+21 might provide another one. 

By the spin angular momentum of the pulsar given in Table 1,3,5, one can easily obtain 
its moment of inertia, Ji = I.5/45, 2.6/45, 0. 44/45 for solution 1,2,3, respectively, (since pulsar 
period is P = 1.5578ms), which arc consistent with the prediction of neutron star equation 
of state. It is interesting that PSR B1937+21 having the smallest pulsar period, P, might 
also has the smallest x of all miUisecond pulsars. 

I would like to thank Prof. T.Lu for constructive suggestions. 1 also thank Prof. Z.R.Wang, 
X.D.Li, Z.Y.Li, Z.D.Dai, D.M.Wei for helpful discussion on the manuscript. 
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Table 1: Solution 1, obtained by fitting the measured data of Fig 1 



mi (Mo) 


m2(MQ) 


e 


n(d) 


i 




1.30 


0.0373 


0.0426 


0.0708 


7.25 • 10-=^ 


0.779 


5*1 (g cm^s-i) 


S2{g cm^s-i) 










1.80 • 10^^ 


1.83 • 10^8 


0.496 


0.0677 


0.983 


2.56 



All angles are in radian. 5*1 corresponds to the moment of inertia of the pulsar, Ii = I.53I45. 



Table 2: Derived parameters through solution 1 



f^o(s-i) 


x{s) 




X 




x(s 1) 








-5.16 ■ 10-1° 


3.73- 


10-^ 


8.33 


■ 10-1^ 


-4.37 


■ 10-22 


3.28- 


10-21 




A 












pGR 




1.11 ■ 10-^ 


1.43 • 


10-11 


2.52 


■ 10-18 


-2.96 


■ 10-27 


8.06- 


10-14 



X and X are calculated by assuming r}Q = ujq 



Table 3: Solution 2 



mi(M0) 


m2(Mo) 


e P(,(d) 


i 


7L(yr) 


1.38 


0.0364 


0.147 0.100 


0.0101 


0.779 


5*1 (g cm2s-i) 


5*2 (g cm2s-i) 




(t>i 




3.04 • 10^8 


1.44 • 10^8 


1.46 -0.602 


3.91 


4.15 



All angles are in radian. Si corresponds to the moment of inertia of the pulsar, Ii = 2.59/45. 



Table 4: Derived parameters through solution 2 





x{s) 




X 




x{s 








-1.94- 10-1° 


6.16- 


10-4 


-6.63- 


10-14 


9.11 


■ 10-22 


2.49 • 


10-21 


c:;«^(s-i) 


Pb 




A(s-i) 


1 




[s-2) 


PGR 




6.59 ■ 10-9 


3.19- 


^ 10-11 


-3.07- 


10-2° 


2.14 


■ 10-28 


5.22 • 


10-14 



X and X are calculated by assuming r}Q = ujq 
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Table 5: Solution 3 



mi(MQ) 


m2(M0) 


e 


A(d) 


i 


To(yr) 


1.88 


0.0126 


0.414 


0.0975 


0.250 


0.779 


5'i(g cm^s"-^) 


S'2(g cm^s"^) 




^JS2 






5.15 • 10^^ 


1.85 • 10 "-^ 


1.31 


0.0171 


3.31 


0.687 



All angles are in radian. Si corresponds to the moment of inertia of the pulsar, /i = O.439/45. 



Table 6: Derived parameters through solution 3 





x{s) 


X 


x(s 1) 


j.GR 




2.95 • 10-1^ 


4.25 • 10-=^ 


-9.90 • 10-16 


1.44 • 10-23 


2.02- 


10-20 


a;G^'(s-i) 


A 


A(s-i) 


Pr(s-2) 






9.86 • 10-*^ 


-4.84 • 10-^2 


3.39 • 10-2° 


-5.62 • 10-29 


5.99- 


10-14 



X and X are calculated by assuming 770 = '^o 




Fig. 1. — The scattered points are measured by Kaspi et al. (1994), the sohd hne is given 
by the orbital precession model. 




Fig. 2. — The scattered points are measured by Kaspi et al. (1994), the sohd hne is given 
by the orbital precession model. 




Fig. 3. — The scattered points are measured by Kaspi et al. (1994), the sohd hne is given 
by the orbital precession modeL 



